package statistics;

import java.math.BigDecimal;
import java.util.Map;
import java.util.concurrent.ConcurrentHashMap;

import pac.KeyGen;

public class ExactTest extends AbstractStatisticTests 
{	
	
	static BigDecimal[] lnFactorialList = new BigDecimal[5000];
	static {
		lnFactorialList[0] = new BigDecimal("0").setScale(100);
		for(int i=1; i<5000; i++) {
			BigDecimal number = new BigDecimal( Math.log( i ) ).setScale(100);
			lnFactorialList[i] = lnFactorialList[i-1].add(number);
		}
	}
	
	private IStatisticalTest asymptoticTest;
	public ExactTest (IStatisticalTest testOfChoice) {
		asymptoticTest = testOfChoice;
	}
	
	private static BigDecimal choose ( int m, int k ) 
	{	
		if (k<0 || k>m)
			return new BigDecimal(1);
		
		BigDecimal numerator = new BigDecimal(1).setScale(100);
		BigDecimal denominator = new BigDecimal(1).setScale(100);
		
		int greater = k >= (m-k) ? k : (m-k); 
		int lesser = k < (m-k) ? k : (m-k); 
		
		for(int i=m; i>greater; i--)
			numerator = numerator.multiply(new BigDecimal(i));
		for(int i=lesser; i>0; i--)
			denominator = denominator.multiply(new BigDecimal(i));
		BigDecimal m_choose_k = numerator.divide(denominator, BigDecimal.ROUND_HALF_UP);
		
		return m_choose_k.setScale(100);
	}
	
	private static int[] configUp (int[] configuration) 
	{
		configuration[0]++;
		configuration[1]--;
		configuration[2]--;
		configuration[3]++;
		return configuration;
	}
	
	private static int[] configDown (int[] configuration) 
	{
		configuration[0]--;
		configuration[1]++;
		configuration[2]++;
		configuration[3]--;
		return configuration;
	}
	
	private int[] minConfig (int[] config) {
		int[] configuration = config.clone();
		int delta = Math.min(configuration[0], configuration[3]);
		configuration[0]-=delta;
		configuration[1]+=delta;
		configuration[2]+=delta;
		configuration[3]-=delta;
		return configuration;
	}
		
	//makes a long key for a combination 
	//e.g.1 : 400c2 is coded as 000400-00000
	//e.g.2 : A 2x2 table with the counts 1900/102/40/5 is coded as 00190-00102-00040-00005
	private static KeyGen keygen = new KeyGen(10000); 	

	
	//See http://www.ibm.com/developerworks/java/library/j-jtp07233/index.html#N10244
	private static volatile Map<Double, BigDecimal> configurationToPointProbability = new ConcurrentHashMap<Double, BigDecimal>();
	private static BigDecimal calculatePointEventProbability (int[] config) 
	{
		int[] configuration = config.clone();
		double key = keygen.makeKey(configuration);
		BigDecimal probability = configurationToPointProbability.get(key);
		
		if ( probability == null ) 
		{
			int m = configuration[0] + configuration[1];
			int N_minus_m = configuration[2] + configuration[3];
			int N = m + N_minus_m;
			int n = configuration[0] + configuration[2];
			int k = configuration[0];
			int n_minus_k = n - k;
			
			//log probability of hyper-geometric distribution
			
			BigDecimal numerator = choose(m, k).multiply( choose(N_minus_m, n_minus_k) );
			BigDecimal denominator = choose(N, n);
			probability = numerator.divide( denominator , BigDecimal.ROUND_HALF_UP ) ;
			
			configurationToPointProbability.put ( key, probability );
		}		
		
		return probability;
	}
	
	
	private static volatile Map<Double, BigDecimal> configurationToCummulativeProbability = new ConcurrentHashMap<Double, BigDecimal>();
	private static BigDecimal calculateCummulativeEventProbability (int[] config) 
	{
		int[] configuration = config.clone();
		double key = keygen.makeKey(configuration);
		BigDecimal cumm_probability = configurationToCummulativeProbability.get(key);
		
		if ( cumm_probability == null ) {
			cumm_probability = new BigDecimal("0");		
			
			while ( isValid(configuration) )  
			{
				BigDecimal probability = calculatePointEventProbability(configuration);
				cumm_probability = cumm_probability.add(probability);
				configuration = configDown(configuration);
			} 
			
			configurationToCummulativeProbability.put(key, cumm_probability);
		}
		
		return cumm_probability;
	}
	
	
	public double TestCasesOnly ( int order, int[][][] observedCases ) 
	{		
		int[] configuration = make1DArray(observedCases);

		BigDecimal pValue = new BigDecimal("0");		
		while( isValid(configuration) ) 
		{
			pValue = pValue.add ( calculatePointEventProbability(configuration) );
			configuration = configUp(configuration);
		}
		
		return pValue.doubleValue();
	}

	
	
	//adds up the events that are more extreme (less likely) than the current case control configuration
	public double TestCasesVsControls ( int order, int[][][] observedCases, int[][][] observedControls ) 
	{		
		
		double boundaryZscore = asymptoticTest.TestCasesVsControls_Zscore ( order, observedCases, observedControls );
		
		// define the side of the distribution you are looking at.
		// if this is control enriched, then swap and test for more extreme control enrichment.
		if (boundaryZscore < 0) {
			int[][][] temp = observedCases;
			observedCases = observedControls;
			observedControls = temp;
			boundaryZscore = - boundaryZscore;
		}
			
		int[] case_configuration = make1DArray(observedCases);
		int[] control_configuration = make1DArray(observedControls);
		
		BigDecimal	exactPvalue = new BigDecimal("0"); //set it to zero, so you can keep count.
		
		double currentZscore = 0.0;
		int[][][] control_config_3D;
		int[][][] case_config_3D;
		
		//increment cases
		case_configuration = minConfig(case_configuration);
		while( isValid(case_configuration) ) 
		{			
			case_config_3D = make3DArray(case_configuration);
			
			control_configuration = minConfig(control_configuration);					
			while( isValid(control_configuration) ) 
			{
				control_config_3D = make3DArray(control_configuration);
				
				currentZscore = asymptoticTest.TestCasesVsControls_Zscore ( order, case_config_3D, control_config_3D );
											
				//you've gone too far. decrement and break.
				if ( currentZscore <  boundaryZscore ) 
					break;	
				
				//while you can get away with it, increment the controls.
				control_configuration = configUp(control_configuration);	
			}
			
			//check if you've gone too far ahead
			control_configuration = configDown(control_configuration);
			if( !isValid(control_configuration) )
				control_configuration = configUp(control_configuration);
			
			//redo some work to see where you stand
			control_config_3D = make3DArray(control_configuration);
			currentZscore = asymptoticTest.TestCasesVsControls_Zscore ( order, case_config_3D, control_config_3D );
			
			//this is as far as you can go by incrementing controls. calculate the exact p
			if ( currentZscore >= boundaryZscore && isValid(case_configuration) && isValid(control_configuration) ) 
			{			
				BigDecimal casePointProbability = calculatePointEventProbability(case_configuration);
				BigDecimal controlCummulativeProbability = calculateCummulativeEventProbability(control_configuration);	

//				Z-scores aren't negative anymore.
//				if ( boundaryZscore < 0 )
//					controlCummulativeProbability = new BigDecimal(1).subtract( controlCummulativeProbability );
				
				BigDecimal configurationPvalue = casePointProbability.multiply( controlCummulativeProbability );
				exactPvalue = exactPvalue.add( configurationPvalue );
			}
			
			//proceed to the next extreme case configuration.
			case_configuration = configUp(case_configuration);
			
		}
		
		return exactPvalue.doubleValue();
	}
	
	
	private static boolean isValid(int[] configuration) 
	{
		if ( configuration[0] >=0 
				&& configuration[1] >=0 
				&& configuration[2] >=0 
				&& configuration[3] >=0 )
			return true;
		else
			return false;
	}
	
}
